Self-consistent calculation of metamaterials with gain 
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We present a computational scheme allowing for a self-consistent treatment of a dispersive metallic photonic 
metamaterial coupled to a gain material incorporated into the nanostructure. The gain is described by a generic 
four-level system. A critical pumping rate exists for compensating the loss of the metamaterial. Nonlinearities 
arise due to gain depletion beyond a certain critical strength of a test field. Transmission, reflection, and ab- 
sorption data as well as the retrieved effective parameters are presented for a lattice of resonant square cylinders 
embedded in layers of gain material and split ring resonators with gain material embedded into the gaps. 
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The field of metamaterials 1 ^ is driven by fascinating 
and far-reaching theoretical visions such as, e.g., per- 
fect lenses, 3 invisibility cloaking, 4,5 and enhanced optical 
nonlinearities. 6 This emerging field has seen spectacular ex- 
perimental progress in recent years. 1 ' 2 Yet, losses are orders of 
magnitude too large for the envisioned applications. Achiev- 
ing such reduction by further design optimization appears to 
be out of reach. Thus, incorporation of active media (gain) 
might come as a cure. The dream would be to simply in- 
ject an electrical current into the active medium, leading to 
gain and hence to compensation of the losses. However, 
experiments on such intricate active nanostructures do need 
guidance by theory via self-consistent calculations (using the 
semi-classical theory of lasing) for realistic gain materials that 
can be incorporated into or close to dispersive media to re- 
duce the losses at THz or optical frequencies. The need for 
self-consistent calculations stems from the fact that increas- 
ing the gain in the metamaterial, the metamaterial properties 
change, in turn changes the coupling to the gain medium until 
a steady-state is reached. A specific geometry to overcome 
the severe loss problem of optical metamaterials and to en- 
able bulk metamaterials with negative magnetic and electric 
response and controllable dispersion at optical frequencies is 
to interleave active optically pumped gain material layers with 
the passive metamaterial lattice. 

For reference, the best fabricated negative-index material 
operating at around 1.4 /im wavelength 7 has shown a figure of 
merit (FOM) = — Re(n) /Im(n) ~ 3, where n is the effective 
refractive index. This experimental result is equivalent to an 
absolute absorption coefficient of a = 3 x 10 4 cm , which 
is even larger than the absorption of typical direct-gap semi- 
conductors such as, e.g., GaAs (where a — 10 4 cm -1 ). So it 
looks difficult to compensate the losses with this simple type 
of analysis, which assumes that the bulk gain coefficient is 
needed. However, the effective gain coefficient, derived from 
self-consistent microscopic calculations, is a more appropri- 
ate measure of the combined system of metamaterial and gain. 
Due to pronounced local-field enhancement effects in the spa- 
tial vicinity of the dispersive metamaterial, the effective gain 
coefficient can be substantially larger than its bulk counter- 



part. While early models using simplified gain-mechanisms 
such as explicitly forcing negative imaginary parts of the local 
gain material's response function produce unrealistic strictly 
linear gain, our self-consistent approach presented below al- 
lows for determining the range of parameters for which one 
can realistically expect linear amplification and linear loss 
compensation in the metamaterial. To fully understand the 
coupled metamaterial-gain system, we have to deal with time- 
dependent wave equations in metamaterial systems by cou- 
pling Maxwell's equations with the rate equations of elec- 
tron populations describing a multi-level gain system in semi- 
classical theory £ 

In this paper, we apply a detailed computational model to 
the problem of metamaterials with gain. The generic four- 
level atomic system tracks fields and occupation numbers at 
each point in space, taking into account energy exchange be- 
tween atoms and fields, electronic pumping and non-radiative 
decays. 8 An external mechanism pumps electrons from the 
ground state level Ao to the third level A3 at a certain pump- 
ing rate r pump , which is proportional to the optical pumping 
intensity in an experiment. After a short lifetime T32 electrons 
transfer non-radiatively into the metastable second level N2 . 
The second level (N2) and the first level (Ni) are called the 
upper and lower lasing levels. Electrons can be transferred 
from the upper to the lower lasing level by spontaneous and 
stimulated emission. At last, electrons transfer quickly and 
non-radiatively from the first level (N{) to the ground state 
level (No). The lifetimes and energies of the upper and lower 
lasing levels are T21, E% and Tin, E\, respectively. The cen- 
ter frequency of the radiation is oj a = (E2 — Ei)/H which is 
chosen to equal 2tt x 10 14 Hz. The parameters T32, t%\, and 
Tin are chosen 5 x 10~ 14 , 5 x 10~ 12 , and 5 x 10~ 14 s, re- 
spectively. The total electron density, No(t = 0) = No(t) + 
Ni(t) + N 2 {t) + N 3 (t) = 5.0 x 10 23 /rn 3 , and the pump 
rate r pump are controlled variables according to the experi- 
ment. The time-dependent Maxwell equations are given by 
V x E = -dB/dtand Vx H = ee Q dE/dt + dP/dt, where 
B = /i/i H and P is the dispersive electric polarization den- 
sity from which the amplification and gain can be obtained. 
Following the single electron case, we can show 8 that the po- 
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FIG. 1: (Color online) Square lattice of dielectric square cylinders 
(blue) that have a Lorentz behavior embedded in layers of gain ma- 
terial (red). The dielectric constant of the cylinders is given by 
e = 1 + ujp/(ujl - 2iLo-y - uj 2 ), where f p = lu p /2iy = 100 THz 
and 7 = 2nf, and / takes different values in the cases we have 
examined. The dimensions are a = 80 nm, wl = 40 nm and 



larization density P(r, t) in the presence of an electric field 
obeys locally the following equation of motion 



2 P(t) , n dP(t) 



dt 2 



r - 



dt 



^P(t) = -a a AN(t)E(t) (1) 



where r a is the linewidth of the atomic transition uj a and is 
equal to 2?r x 5 x 10 12 Hz or 2tt x 20 x 10 12 Hz. The factor 
AN(r,t) — N2(r,t) — N%(r,t) is the population inversion 
that drives the polarization, and a a is the coupling strength 
of P to the external electric field and its value is taken to be 
10~ 4 C 2 /kg. It follows 8 from Eqn. Q]that the amplification 
line shape is Lorentzian and homogeneously broadened.— The 
occupation numbers at each spatial point vary according to 
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where 



huj, 



■E 



^ is the induced radiation rate or excitation 



rate depending on its sign. 

In order to solve the behavior of the active materials in the 
electromagnetic fields numerically, the finite-difference time- 
domain (FDTD) technique is utilized, 10 using an approach 
similar to the one outlined in Refs. 10-12. In the FDTD cal- 
culations, the discrete time and space steps are chosen to be 
At = 8.33 x 10~ 18 s and Ax = 5.0 x 10~ 9 m for simulations 
on the structure as shown in Fig. 1, and At = 8.33 x 10" 19 s 
and Ax — 1.0 x 10 _9 m for simulations on the structure as 
shown in Fig. 5. The initial condition is that all the electrons 
are in the ground state, so there is no field, no polarization 
and no spontaneous emission. Then the electrons are pumped 
from No to N3 (then relaxing to A^) with a constant pump- 
ing rate r pump . The system begins to evolve according to the 
system of equations above. 

We have performed numerical simulations on one- 
dimensional (ID) and two-dimensional (2D) systems with 
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FIG. 2: (Color online) The transmittance vs. probe field amplitude 
for the loss-compensated metamaterial of Fig. 1 with gain bandwidth 
5 THz, loss bandwidth 20 THz (i.e., / = 10 THz), for different 
pumping rates r pump . r pump is increased from 2.15 x 10 9 s _1 
(lowest) to 3.05 x lO'-'s" 1 (highest) in steps of 0.1 x 10 9 s~\ 
The metamaterial response is linear in a very wide range. When 
the loss-compensated transmission is exactly unity, the pumping rate 
r pump = 2.65 x 10 9 s _1 , which is called the critical pumping rate. 
For incident fields stronger than 10 4 V/m this metamaterial becomes 
non-linear. 



gain^ Previous studie s 14 ' 15 ' 16 ' 17 ' 18 have considered loss re- 
duction by incorporating gain but where not self-consistent 
(see introduction) . 14 ' 15 ' 16 ' 17 As the first simple model system, 
we will discuss a 2D metamaterial system (shown in Fig. 1) 
which consists of layers of gain material and dielectric wires 
that have a resonant Lorentz type electric response to emu- 
late the resonant elements in a realistic metamaterial. We will 
have to study whether we will be able to compensate the losses 
of the metamaterials associated with the Lorentz resonance in 
the wires by the amplification provided by the gain material 
layers without destroying the linear response of the metama- 
terial. First we generate a narrow band Gaussian pulse of a 
given amplitude and let it propagate through the metamaterial 
without gain, and we calculate the transmitted signal emerging 
from the metamaterial which has also Gaussian profile but the 
amplitude is much smaller than that of the incident pulse de- 
pending on how much loss occurs in the metamaterial. Then 
we introduce the gain and start increasing the pumping rate 
and find a critical pumping rate, r pump = 2.65 x 10 9 s _1 , for 
which the transmitted pulse is of the same amplitude as the 
incident pulse. In addition, for fixed pumping rate, we start 
increasing the amplitude of the incident Gaussian pulse and 
we would like to see how high we can go in the strength of 
the incident electric field and still have full compensation of 
the losses, i.e. the transmitted signal equals the incident sig- 
nal, independent on the signal strength. In this region we have 
compensated loss and still linear response of the metamaterial; 
here, the shape of the transmitted Gaussian is only affected by 
the dispersion but not dependent on the signal strength. 

We have calculated the transmission versus the strength of 
the electric field of the incident signal for several pumping 
rates close to the critical pumping rate. As shown in Fig. 2, we 
found that for a rather broad region of low intensity input sig- 
nal we have a linear response all the way up to incident elec- 
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FIG. 3: (Color online) The critical pumping rates for different num- 
bers of layers of the system in Fig. 1 . Parameters for gain and dielec- 
tric materials are the same as Fig. 2. 
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FIG. 5: (Color online) Geometry for a unit cell of the square SRR 
system with gain. The gain (shown in orange) is introduced in the 
gap region of the SRR. The dimensions are a — 100 nm, I = 80 nm, 
t = 5 nm, d — 4 nm and w = 15 nm. 
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FIG. 4: (Color online) The retrieved results for the real and the imag- 
inary parts of the effective permittivity, e, with and without gain. In 
addition, we have plotted Im(e 9 ) versus frequency. Below compen- 
sation, t — 0.75; Gain and Lorentz bandwidths are 20 and 5 THz, 
respectively. 



trie field of 10 3 V/m. If we use only three layers (rods - gain 
material - rods), the critical pumping rate is 4.85 x 10 9 s _1 , 
which is two times higher than the 19-layer case of Fig. 1. 
In Fig. 3, we present detailed results for the critical pump- 
ing rate versus the number of layers of the system shown in 
Fig. 1. Notice that as the number of layers increases, the crit- 
ical r pump decreases. The linear regime for three layers ex- 
ists up to 10 4 V/m, and for higher strength drops slower than 
that of Fig. 2. In all the following simulations, the strength 
of the incident signal is chosen to be 10 V/m, which is far 
away from 10 3 V/m, so we operate in the linear regime of the 
metamaterial. As an example, we have studied three layers, 
rods - gain material - rods, to see how much r pump we need 
to compensate the losses. As expected, we found that r pump 
is proportional to the imaginary part of the permittivity e of 
the dielectric. 

We first present results for three layers of the system shown 
in Fig. 1. First, the full width at half maximum (FWHM) for 
Lorentz dielectric and gain are chosen to be 5 and 20 THz, re- 
spectively. With the introduction of gain the absorption at the 
resonance frequency of 100 THz decreases, ultimately reach- 



ing zero (not shown). So the gain compensates the losses. 
In Fig. 4, we plot the retrieved results for the real and imagi- 
nary parts of e without gain and with gain slightly below com- 
pensation (see Ref. 19 for the retrieval method). Notice that 
we can have the Re(e) rj -1 with Im(e) w at 102 THz, 
slightly off the resonance frequency. From Fig. 4, one can 
also see that Re(e) » 2.5 with Im(e) w at 97 THz. So 
one can obtain a lossless metamaterial with positive or neg- 
ative Re(e). Once we introduce gain, the imaginary part of 
e of our total system with gain is equal to the sum of Im(e) 
without gain and the imaginary part of e g , the dielectric func- 
tion of the gain material. This result is unexpected, because 
there is no coupling between the 2D Lorentz dielectric with 
the gain material. This is indeed true because of the continu- 
ous shape of the Lorentz dielectric cylinders and the gain ma- 
terial slabs have zero depolarization field. In contrast to finite 
length wires (hence a 3D problem) where the dipole interac- 
tions between Lorentz dielectric and gain material would be 
dominated by the quasi-static nearfield 0(1 /r 3 ), here the in- 
teraction is order 0(u> In \ kr\), only via the propagating field, 
and much weaker. Therefore, for this 2D model, gain and 
loss are approximately independent. The behavior would ob- 
viously be different in a 3D situation, which, however, is com- 
putationally excessively demanding. Thus, we consider a 2D 
version of the split ring resonator (SRR) as a more realistic 
and also more relevant model. Here, the relevant polarization 
is across the finite SRR gap and, therefore, the coupling to the 
gain material is in fact dipole like. 

In Fig. 5, we present the unit cell of our SRR system with 
gain material embedded in the SRR gap. The dimensions of 
the SRR are chosen such that a magnetic resonance frequency 
at 100 THz results, which can overlap with the peak of the 
emission of the gain material. The FWHM of the gain ma- 
terial is 20 THz, and r pump is 1.4 x 10 9 s _1 . Simulations 
are done for one layer of the square SRR. In Fig. 6a, we plot 
the retrieved results of the real and the imaginary parts of the 
magnetic permeability, /x, with and without gain. With the in- 
troduction of gain, the weak and broad resonant effective fi 
(FWHM = 5.85 THz) of the lossy SRR becomes strong and 
narrow (FWHM = 1.66 THz); the gain effectively undamps 
the LCR resonance of the SRR. Notice that here losses in the 
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FIG. 6: (Color online) The retrieved results for the real and the imag- 
inary parts of (a) the effective permeability fi and (b) the correspond- 
ing effective index of refraction n, with and without gain for a pump- 
ing rate r pump = 1.4 x 10 9 s" 1 . The gain bandwidth is 20 THz. 
Notice that the width of the resonance with gain is 1.66 THz. 

magnetic effective response are compensated by electric gain 
in the SRR gap. So with the introduction of gain, we obtain 
a negative fi with a very small imaginary part in an otherwise 
typical SRR response, which means that the losses have been 
compensated by the gain. In Fig. 6b, we plot the retrieved re- 
sults for the effective index of refraction n, with and without 
gain. Note that for a lossless SRR n is purely real away from 
the resonance and imaginary in a small band above the reso- 
nance where fi is negative. Comparing Re(n) slightly below 
the resonance at 97 THz, we find an effective extinction co- 
efficient a = (w/c) Im(n) 3.50 x 10 4 cm _1 without gain, 
and a 1.24 x 10 4 cm _1 with gain, and hence an effective 
amplification of a w —2.26 x 10 4 cm -1 . This is much larger 



than the expected amplification a « — 1.39 x 10 3 cm _1 for 
the gain material at the given pumping rate. 20 The difference 
can be explained by the field enhancement in the gap of the 
resonant SRR. The induced electric field in the gap is around 
550 V/m, which is still in the linear regime, and the incident 
electric field is 10 V/m. Indeed, taking the observed field en- 
hancement factor in the SRR gap of w 55, the energy per unit 
cell produced by the gain material inside the gap is w 18 times 
larger than for the homogeneous gain medium which com- 
pares very well to the factor « 20 between the simulated SRR 
effective medium and the homogeneous gain medium. If we 
further increase the pumping rate the magnetic resonance be- 
comes even narrower (0.96 THz for r pump = 1.8 x 10 9 s _1 ). 
When the pumping rate reaches r pump = 1.9 x 10 9 s _1 , 
Im(/Lt) becomes negative and we have overcompensated at 
the resonance frequency. By increasing T pump even more 
(«5x 10 9 s _1 ) one starts seeing lasing (spasing) 21 ' 22 in our 
system (not shown), which is not in the focus of this work. As 
long as we are in the linear regime, we do not need to have 
a self-consistent calculation, our results agree very well with 
the results obtained using the susceptibilities given in Ref. 9. 
However, the self-consistent calculation is necessary to deter- 
mine the range of signals for which we can expect approxi- 
mately linear response and it is needed if we have very strong 
fields and we are in the nonlinear regime, especially when we 
want to study lasing. 

In conclusion, we have proposed and numerically solved 
a self-consistent model incorporating gain in 2D dispersive 
metamaterials. We show numerically that one can compen- 
sate the losses of the dispersive metamaterials. There is a rel- 
atively wide range of signal amplitudes for which the loss- 
compensated metamaterial still behaves linearly; at higher 
amplitudes the response is non-linear due to the gain. As an 
example, we have demonstrated that the losses of the mag- 
netic susceptibility \i of the SRR can be easily compensated 
by the gain material. The pumping rate needed to compensate 
the loss is much smaller than the bulk gain. This aspect is due 
to the strong local-field enhancement inside the SRR gap. 
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